home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMATNL.H < prev    next >
C/C++ Source or Header  |  1995-01-11  |  8KB  |  232 lines
  1. //$$ newmatnl.h           definition file for non-linear optimisation
  2.  
  3. // Copyright (C) 1993,4,5: R B Davies
  4.  
  5. #ifndef NEWMATNL_LIB
  6. #define NEWMATNL_LIB 0
  7.  
  8. #include "newmat.h"
  9.  
  10. /*
  11.  
  12. This is a beginning of a series of classes for non-linear optimisation.
  13.  
  14. At present there are two classes. FindMaximum2 is the basic optimisation
  15. strategy when one is doing an optimisation where one has first
  16. derivatives and estimates of the second derivatives. Class
  17. NonLinearLeastSquares is derived from FindMaximum2. This provides the
  18. functions that calculate function values and derivatives.
  19.  
  20.  
  21.  
  22.    class FindMaximum2
  23.  
  24. Suppose T is the ColumnVector of parameters, F(T) the function we want
  25. to maximise, D(T) the ColumnVector of derivatives of F with respect to
  26. T, and S(T) the matrix of second derivatives.
  27.  
  28. Then the basic iteration is given a value of T, update it to
  29.  
  30.            T - S.i() * D
  31.  
  32. where .i() denotes inverse.
  33.  
  34. If F was quadratic this would give exactly the right answer (except it
  35. might get a minimum rather than a maximum). Since F is not usually
  36. quadratic, the simple procedure would be to recalculate S and D with the
  37. new value of T and keep iterating until the process converges. This is
  38. known as the method of conjugate gradients.
  39.  
  40. In practice, this method may not converge. FindMaximum2 considers an
  41. iteration of the form
  42.  
  43.            T - x * S.i() * D
  44.  
  45. where x is a number. It tries x = 1 and uses the values of F and its
  46. slope with respect to x at x = 0 and x = 1 to fit a cubic in x. It then
  47. choses x to maximise the resulting function. This gives our new value of
  48. T. The program checks that the value of F is getting better and carries
  49. out a variety of strategies if it isn't.
  50.  
  51. The program also has a second strategy. If the successive values of T
  52. seem to be lying along a curve - eg we are following along a curved
  53. ridge, the program will try to fit this ridge and project along it. This
  54. doesn't work at present and is commented out.
  55.  
  56. FindMaximum2 has three virtual functions which need to be over-ridden by
  57. a derived class.
  58.  
  59.    void Value(const ColumnVector& T, Boolean wg, Real& f, Boolean& oorg);
  60.  
  61. T is the column vector of parameters. The function returns the value of
  62. the function to f, but may instead set oorg to TRUE if the parameter
  63. values are not valid. If wg is TRUE it may also calculate and store the
  64. second derivative information.
  65.  
  66.    Boolean NextPoint(ColumnVector& H, Real& d);
  67.  
  68. Using the value of T provided in the previous call of Value, find the
  69. conjugate gradients adjustment to T, that is - S.i() * D. Also return
  70.  
  71.            d = D.t() * S.i() * D.
  72.  
  73. NextPoint should return TRUE if it considers that the process has
  74. converged (d very small) and FALSE otherwise. The previous call of Value
  75. will have set wg to TRUE, so that S will be available.
  76.  
  77.    Real LastDerivative(const ColumnVector& H);
  78.  
  79. Return the scalar product of H and the vector of derivatives at the last
  80. value of T.
  81.  
  82. The function Fit is the function that calls the iteration.
  83.  
  84.    void Fit(ColumnVector&, int);
  85.  
  86. The arguments are the trial parameter values as a ColumnVector and the
  87. maximum number of iterations. The program calls a DataException if the
  88. initial parameters are not valid and a ConvergenceException if the
  89. process fails to converge.
  90.  
  91.  
  92.    class NonLinearLeastSquares
  93.  
  94. This class is derived from FindMaximum2 and carries out a non-linear
  95. least squares fit. It uses a QR decomposition to carry out the
  96. operations required by FindMaximum2.
  97.  
  98. A prototype class R1_Col_I_D is provided. The user needs to derive a
  99. class from this which includes functions the predicted value of each
  100. observation its derivatives. An object from this class has to be
  101. provided to class NonLinearLeastSquares.
  102.  
  103. Suppose we observe n normal random variables with the same unknown
  104. variance and such the i-th one has expected value given by f(i,P)
  105. where P is a column vector of unknown parameters and f is a known
  106. function. We wish to estimate P.
  107.  
  108. First derive a class from R1_Col_I_D and override Real operator()(int i)
  109. to give the value of the function f in terms of i and the ColumnVector
  110. para defined in class R1_CoL_I_D. Also override ReturnMatrix
  111. Derivatives() to give the derivates of f at para and the value of i
  112. used in the preceeding call to operator(). Return the result as a
  113. RowVector. Construct an object from this class. Suppose in what follows
  114. it is called pred.
  115.  
  116. Now constuct a NonLinearLeastSquaresObject accessing pred and optionally
  117. an iteration limit and an accuracy critierion.
  118.  
  119.    NonLinearLeastSquares NLLS(pred, 1000, 0.0001);
  120.  
  121. The accuracy critierion should be somewhat less than one and 0.0001 is
  122. about the smallest sensible value.
  123.  
  124. Define a ColumnVector P containing a guess at the value of the unknown
  125. parameter, and a ColumnVector Y containing the unknown data. Call
  126.  
  127.    NLLS.Fit(Y,P);
  128.  
  129. If the process converges, P will contain the estimates of the unknow
  130. paramters. If it doesn't converge an exception will be generated.
  131.  
  132. The following member functions can be called after you have done a fit.
  133.  
  134. Real ResidualVariance() const;
  135.  
  136. The estimate of the variance of the observations.
  137.  
  138. void GetResiduals(ColumnVector& Z) const;
  139.  
  140. The residuals of the individual observations.
  141.  
  142. void GetStandardErrors(ColumnVector&);
  143.  
  144. The standard errors of the observations.
  145.  
  146. void GetCorrelations(SymmetricMatrix&);
  147.  
  148. The correlations of the observations.
  149.  
  150. void GetHatDiagonal(DiagonalMatrix&) const;
  151.  
  152. Forms a diagonal matrix of values between 0 and 1. If the i-th value is
  153. larger than, say 0.2, then the i-th data value could have an undue
  154. influence on your estimates.
  155.  
  156.  
  157. */
  158.  
  159. class FindMaximum2
  160. {
  161.    virtual void Value(const ColumnVector&, Boolean, Real&, Boolean&) = 0;
  162.    virtual Boolean NextPoint(ColumnVector&, Real&) = 0;
  163.    virtual Real LastDerivative(const ColumnVector&) = 0;
  164. public:
  165.    void Fit(ColumnVector&, int);
  166. };
  167.  
  168. class R1_Col_I_D
  169. {
  170.    // The prototype for a Real function of a ColumnVector and an
  171.    // integer.
  172.    // You need to derive your function from this one and put in your
  173.    // function for operator() and Derivatives() at least.
  174.    // You may also want to set up a constructor to enter in additional
  175.    // parameter values (that won't vary during the solve).
  176.  
  177. protected:
  178.    ColumnVector para;                 // Current x value
  179.  
  180. public:
  181.    virtual Boolean IsValid() { return TRUE; }
  182.                                        // is the current x value OK
  183.    virtual Real operator()(int i) = 0; // i-th function value at current para
  184.    virtual void Set(const ColumnVector& X) { para = X; }
  185.                                        // set current para
  186.    Boolean IsValid(const ColumnVector& X)
  187.       { Set(X); return IsValid(); }
  188.                                        // set para, check OK
  189.    Real operator()(int i, const ColumnVector& X)
  190.       { Set(X); return operator()(i); }
  191.                                        // set para, return value
  192.    virtual ReturnMatrix Derivatives() = 0;
  193.                                        // return derivatives as RowVector
  194. };
  195.  
  196.  
  197. class NonLinearLeastSquares : public FindMaximum2
  198. {
  199.    // these replace the corresponding functions in FindMaximum2
  200.    void Value(const ColumnVector&, Boolean, Real&, Boolean&);
  201.    Boolean NextPoint(ColumnVector&, Real&);
  202.    Real LastDerivative(const ColumnVector&);
  203.  
  204.    Matrix X;                         // the things we need to do the
  205.    ColumnVector Y;                   // QR triangularisation
  206.    UpperTriangularMatrix U;          // see the write-up in newmata.txt
  207.    ColumnVector M;
  208.    Real errorvar, criterion;
  209.    int n_obs, n_param;
  210.    const ColumnVector* DataPointer;
  211.    RowVector Derivs;
  212.    SymmetricMatrix Covariance;
  213.    DiagonalMatrix SE;
  214.    R1_Col_I_D& Pred;                 // Reference to predictor object
  215.    int Lim;                          // maximum number of iterations
  216.  
  217. public:
  218.    NonLinearLeastSquares(R1_Col_I_D& pred, int lim=1000, Real crit=0.0001)
  219.       : Pred(pred), Lim(lim), criterion(crit) {}
  220.    void Fit(const ColumnVector&, ColumnVector&);
  221.    Real ResidualVariance() const { return errorvar; }
  222.    void GetResiduals(ColumnVector& Z) const { Z = Y; }
  223.    void GetStandardErrors(ColumnVector&);
  224.    void GetCorrelations(SymmetricMatrix&);
  225.    void GetHatDiagonal(DiagonalMatrix&) const;
  226.  
  227. private:
  228.    void MakeCovariance();
  229. };
  230.  
  231. #endif
  232.